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We present a Markov Chain Monte Carlo global analysis of neutrino parameters using both cos- 
mological and experimental data. Results are presented for the combination of all presently available 
data from oscillation experiments, cosmology, and neutrinoless double beta decay. In addition we 
explicitly study the interplay between cosmological, tritium decay and neutrinoless double beta de- 
cay data in determining the neutrino mass parameters. We furthermore discuss how the inference 
of non-neutrino cosmological parameters can benefit from future neutrino mass experiments such as 
the KATRIN tritium decay experiment or neutrinoless double beta decay experiments. 
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I. INTRODUCTION 

The question of neutrino mass is one of the most pro- 
found in modern particle physics. Most plausible models 
of neutrino mass solve the puzzle of why neutrino masses 
are so small by introducing a new scale at high energy, 
and precision studies of neutrino physics therefore hold 
the potential to investigate physics at scales beyond those 
reachable in current accelerator experiments. They also 
make the study of the possible Majorana nature of neu- 
trinos possible (see [l|, [U, 9 for a thorough discussion 
of this). While the neutrino mass differences have now 
been measured at about 10% precision by oscillation ex- 
periments (see e.g. 0, [H) the absolute mass scale remains 
unknown and inaccessible to oscillation experiments. 

There are, however, several possible paths to measur- 
ing the absolute neutrino mass. The kinematical effect 
of neutrino mass can be probed either via its effect on 
the beta decay spectrum or via its effect on cosmological 
structure formation. If neutrinos are Majorana particles 
a different possibility is to search for neutrinoless dou- 
ble beta decay because the transition probability for this 
process is proportional to the neutrino mass squared. 

In the past year there have been several papers dis- 
cussing how to unify the data analysis for the various 
approaches @, @- This is a no n- trivial issue, given that 
completely different physics is involved and that the three 
probes are actually sensitive to three distinct observablcs. 

Here we present a new Markov Chain Monte Carlo 
global analysis of neutrino parameters using both cos- 
mological and experimental data. The analysis soft- 
ware is based on the CosmoMC Markov Chain Monte 
Carlo (MCMC) package for cosmological parameter esti- 
mation 0, Q , appropriately modified to incorporate all 
parameters related to neutrino physics. This approach 
uses Baycsian inference instead of the frcquentist method 
commonly used in particle physics. The approach is 
somewhat similar to the MCMC technique developed in 
Q to constrain MSSM parameters. However, a key differ- 
ence is that here we keep the full cosmological parameter 
estimation which allows for a closer study of the inter- 
play between neutrino data and cosmological parameter 



estimation. 

In Section II we describe the methodology used and 
in Section III we present the main results for various 
different assumptions about present and future data, as 
well as different parameter spaces. Finally we present 
our conclusions in Section IV. 



II. METHODOLOGY 

The MCMC Bayesian inference approach has been de- 
scribed in detail for instance in 0, U3] ■ Based on assumed 
priors on each parameter it samples the likelihood func- 
tion using the Markov Chain Monte Carlo method and 
from that the posterior credible intervals for all parame- 
ters can be calculated. Before running the Markov chains 
it is therefore necessary to specify both the parameters 
to be used and the priors on all parameters. 

For the neutrino physics part we have used the mass 
of the lightest eigenstate mj, (mj, = m\ for the normal 
hierachy, mj, = for the inverted hierarchy), the two 
mass differences Amf 2 , A772I3 and the three mixing an- 
gles i?i2, f?23, $13- We also assume that neutrinos are 
Majorana particles so that there are two additional Ma- 
jorana phases <p2, <j>3 \ which together with the mass 
differences and mixing angles specify the observables re- 
lated to absolute neutrino mass (assuming only active 
neutrinos). In total there are then 8 parameters related 
to the neutrino sector. 

There are three separate observables related to the 
three different types of probes. Cosmology is only sensi- 
tive to the neutrino mass, and until the accuracy reaches 
the 0.05 eV level only to the sum of neutrino masses (see 
[Til . IT2I for a thorough discussion of this) 

m v = mi+m2 + ra-3. (1) 



1 These parameters are only important when neutrinoless double 
beta decay data is used 
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In terms of the parameters used in CosmoMC this cor- 
responds to Wi, Am^ and Am^. These are therefore 
the only parameters which are regarded as "slow" in the 
sense that a change in one of them requires recalculation 
of the transfer function for cosmological perturbations. 

At the projected level of accuracy of KATRIN the 
change in the electron energy spectrum can be described 
using a single effective mass parameter which is essen- 
tially the incoherent sum (see e.g. [l3[) 

m f3 = (c?3 c ?2 TO ? + c 2 13 s 2 12 m 2 2 + s 2 3 m 2 ) 1/2 . (2) 

The parameter actually measured in such experiments is 
in fact rrrg which, being a fit-parameter, can be positive 
or negative when measured. 

Conversely, the effective mass measured in neutrinolcss 
double beta decay is the coherent sum [l4|, [3 

m W = \c\ 3 c\ 2 mi + c 2 3 s 2 2 m 2 e l<fe + s 2 3 m 3 e l<fe 1 , (3) 

which allows for phase cancelation. 

The actual parameter measured in any neutrinoless 
double beta decay experiment is the half-life T\/ 2 which 
is related tom^ via the relation [l|| 

J_ = G^|M°f m%, (4) 

1 1/2 

where G 0l/ is a phase-space factor and |M 0!y | 2 the nu- 
clear matrix element squared. In principle the MCMC 
code should use the measured Ti/ 2 and the calculated 
matrix element (both including uncertainties) as param- 
eters instead of m^. However, for simplicity we assume 
a Gaussian error on vnnn with the estimated error on the 
matrix element from [li| [l?], EH ■ 

In the following we have used cosmological parameters 
consistent with the "Vanilla" ACDM model: fi b /i 2 , the 
physical baryon density, £l c h 2 , the physical CDM den- 
sity, A s , the amplitude of primordial fluctuations, n s the 
scalar spectral index, r, the optical depth to reionization, 
and Hq, the Hubble parameter. Spatial flatness has been 
assumed so that Da = 1 — Sl c — Of, — f2„. In total there 
are then 6 parameters related to the cosmological model. 

For the purpose of parameter estimation the publicly 
available CosmoMC Markov chain monte carlo has been 
modified to perform parameter estimation in this 14- 
dimensional parameter space. CosmoMC has been set 
to use the fast /slow parameter scheme. In Table fl] we 
give the list of parameters as well as their priors. Note 
that we also include the two parameters w, the dark en- 
ergy equation of state, and a s , the running of the scalar 
spectral index, both of which are discussed in Section 

mm] 

For the observables related to neutrino oscillation data 
we make the simple assumption of Gaussian errors, given 
by the combination of different experiments. We note 
that this assumption can easily be changed in the code 
and replaced with the full likelihood calculation from ex- 
perimental data. 
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(4- 


12) x 1(T 5 eV 2 
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slow 
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(1- 
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slow 
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Top Hat 
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Top Hat 
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Top Hat 


fast 
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n b h 2 
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Top Hat 
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Top Hat 
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Top Hat 
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Top Hat 
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log(10 10 yl s ) 
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Top Hat 
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0.3-1 


Top Hat 


slow 


w* 




-2 -0 


Top Hat 


slow 






-0.2 - 0.2 


Top Hat 


fast 



TABLE I: Parameters and priors used in the likelihood anal- 
ysis. Cosmological parameters marked with * are used only 
in Section IIITdI 

We use the same constraints as in Q , given by 

Am§ x = (7.92 ± 0.71) x 10 _5 eV 2 

Am 32 = (2.6T.Q jg) x 10~ 3 eV 2 

sin 2 V\ 2 = 0.314±8:gg (5) 

sin 2 tf 2 3 = 0.45±8;5« 

sin 2 tf 2 3 < 0.03, 

with all errors being 2a. The assumption of Gaussian 
errors does not significantly alter any of our results. 



III. RESULTS 

Based on the approach described in the previous sec- 
tion we have calculated the present bound on neutrino 
properties using various combinations of data sets from 
cosmology, tritium decay and neutrinoless double beta 
decay respectively. 



A. Cosmological data 

Cosmological constraints on ^ m„ have been calcu- 
lated by many different authors for various assumptions 
abo ut p arameters and using different data sets (see e.g. 

0, am, sin, nim). 

Here we present just one particular example, which is 
exactly the same as used in (2(|. We use the WMAP 
CMB temperature and polarisation data [27], HH, IpSll . the 
SDSS-LRG and 2dF large scale structure data |30L l3ll . 
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HH , the SDSS-LRG baryon acoustic oscillation data [H| , 
and the SNI-a data set compiled in I34|. Details about 
the cosmological data can be found in [26| . 

Using only the cosmological data we find a bound of 

fn u < 0.50 eV for the minimal ACDM model. We note 
that this is slightly lower than the 0.6 eV found in [26[ for 
the same model and data. This is worth noticing because 
the prior on mj, is logarithmic in the present study while 
it was linear in [2(|. A logarithmic prior on tul tends to 
favour small mj, values because of the large parameter 
space volume at negative log 10 mi and therefore shifts 
the allowed region slightly down. 

This phenomenon is an integral part of Bayesian infer- 
ence because a prior probability distribution needs to be 
specified. In frequentist statistics this problem does not 
occur and the result does not depend on any priors. It 
should be noted that in the limit of Gaussian statistics 
the two methods yield exactly the same result. 

The phenomenon has been recently been studied in the 
context of neutrino properties. For example it was shown 
in (3sT | that Bayesian inference and likelihood maximisa- 
tion give very different results for cosmological parame- 
ters such as the radiation density as long as the likeli- 
hood function is non-Gaussian. As more data is added 
and the likelihood function approaches a Gaussian the 
two methods converge. The question of Bayesian ver- 
sus frequentist statistics was studied in [|| in the context 
of KATRIN. For example the difference between a lin- 
ear and a logarithmic prior on mp was investigated and 
found to have some (not crucial) effect. In conclusion, 
assumptions about priors will have an effect on the pos- 
terior distributions as long as the likelihood function is 
non-Gaussian which is the case for parameters which are 
not extremely well constrained. 

B. Neutrinoless double beta decay 

The upper bound on the effective neutrino mass pro- 
vided by the Heidelberg-Moscow (HM) experiment pro- 
vides an additional and comparable constraint on the ab- 
solute neutrino mass scale [36| . 

We use the constraint 

mpp < 0.27 eV (90%), (6) 

based on the nuclear matrix element calculation in [l7j ]. 
Note that this mass range is more restrictive than what 
was used in [j| because the theoretically predicted half- 
life has been corrected downwards in [l7J compared to 
fl6| . We stress again that the conversion of half- life to 
effective mass 171/3/3 depends strongly on the nuclear ma- 
trix element and that the bound used here could turn out 
to be too restrictive. 

As can be seen from Table UTT adding the HM data does 
shift the allowed range on ^ m„ and m/3/3 down. Since 
the best fit cosmological model in any case has ^ m v = 
it has no influence on other parameters such as fl c h 2 and 
n s . 



Note that we have not derived any cosmological con- 
straint based on the claimed positive evidence from 
Heidelberg-Moscow [37], [H, H^] . Using the same assump- 
tions as above on the nuclear matrix element the claimed 
evidence translates roughly into 0.25 eV < mpp < 0.5 eV 
(90% C.L.). 

C. Future constraints - KATRIN and GERDA 

To get a better idea about the future interplay between 
the three different methods for measuring the absolute 
mass scale we have performed similar likelihood analy- 
ses for the presently available cosmological data together 
with forecasts for the KATRIN beta decay experiment 
[40L |4~0 ] and the GERDA neutrinoless double beta decay 
experiment [H, |43j]. 

For KATRIN we assume a Gaussian la error on to^ 

of ct(to 2 ) = 0.025 eV 2 , roughly in accordance with what 
was used in 0, [44[ . 

For the GERDA neutrinoless double beta decay ex- 
periment we assume a Gaussian error on TTia^ of 0.01 
eV 2 , corresponding roughly to GERDA phase 2 [42j . We 
note that other neutrinoless double beta d ecay experi- 
ments such as MAJORANA and CUORE [II HI will 
reach roughly the same sensitivity in a broadly compa- 
rable time- frame (see e.g. [H[). 

Using these very rough experimental characteristics of 
the two experiments we have proceeded to calculate con- 
straints on neutrino parameters using two different as- 
sumptions: 

(a) In this case we assume no positive detection from 
cither experiment so that the best fit values are mp = 

111/3/3 = 0, 

(b) Here we assume a positive detection from both ex- 
periments, m 2 , = 0.079 eV 2 and m|g = 0.032 eV 2 . 

The last case would for instance be realised in a model 
with normal hierarchy, mj, = nil = x. For both cases 
we perform parameter estimation assuming both normal 
and inverted hierarchy. 

In the first case (a) we note that KATRIN alone would 
significantly tighten the cosmological constraint on ^2 m v 
and GERDA would improve this even further. 171/3/3 is, 
as expected, mainly constrained by adding the GERDA 
data. 

The second case (b) is more interesting from the per- 
spective of combining data sets. The best fit values 
both correspond to roughly 2a evidence for non-zero mp 
and TO/3/3 respectively. However, with the combined data 
771/3/3 = is excluded at roughly 4.5er, likewise mj, = is 
excluded at a similar significance. 

This exercise clearly shows the advantage of analysing 
all neutrino parameters in this global way, instead of sim- 
ply adding constraints. It should be noted that since 
Y^, 7n v ~ is the best fit to present cosmological data the 
case (b) has a best-fit \ 2 which is higher than case (a) 
by A\ 2 = 3.1, i.e. it gives a slightly (not substantially) 
worse fit to cosmological data. 
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Parameter 


Cosmo 


Cosmo+HM 


Cosmo+KATRIN 


Cosmo+KATRIN+GERDA 


Normal hierarchy 






7710,0 = 


«i00,o = 


log 10 rriL (eV) 


1.301 -2.00 


0.836 
L.OIO _ 2 .00 


-1 376 ~°' 892 

l.O/O -2.00 


-1.438 zS:Sg a 


E»» (eV) 


n 91fi 0.489 
U.ZIO 0.080 


0.205 o!o80 


0.189 8:o802 


0-168 0:0790 




o.27i 8;!1J 


o.27i 8i?i 


0.268 8:?72 


0.260 8:?il 


77700 (eV) 


0.0460 o^o 8 


0.0414 Wo 


0.0392 8:o8 x 


0.0317 8'8o 98 


a, ft 2 


o.m giSI 


0.1118:1s 


o.m 8:16:1 


o.m 8:iSI 


n s 


0.951 8;§g 


o.95i8:il 


0.951 8:88 


0.951 8:8i 


Inverted hierarchy 






771/3,0 = 


"i00,o = 


log 10 niL (eV) 


1 or -0.792 
l.OO _ 2 .oo 


i on -0.856 
!.Oi> _2.00 


-1.41 I^o 1 


1 A7 -0.980 
-2.00 


E™» (eV) 


0.232 Ifil 


o.2i28;lf 6 


0.199 8:38 


0.182 ° itl 




A QQi 0.446 
U.OOl o.338 


0.385 8:340 


0.386 8:tS8 


0.394 8:111 


mpp (eV) 


0-0611 g'.oBlo 


0.0554 8: 02 i 8 


0.0538 8:020 


0.0477 8:019 




0.1118:1s; 


o.m 8ioI 


o.m 8:1SI 


o.m 8:1JI 


n 3 


0.951 r 9 !? 


0.951 mi 


o.95i 8:i 2 


o.95i 8:11 


Normal hierarchy 






7710,0 = 0.28 eV, 77100,0 = 0.18 eV 


log 10 m L (eV) 


X 


X 


X 


-0.660 Z° f l 


£m„ (eV) 


X 


X 


X 


0.6748111 




X 


X 


X 


0.330 8:il 


mpp (eV) 


X 


X 


X 


0.166 Ht^ 


Q c 7i 2 


X 


X 


X 


n iic 0.122 
U.110 0.108 


n s 


X 


X 


X 


0.950 8:1!? 


Inverted hierarchy 






77i0,o = 0.28 eV, 77100,0 = 0.18 eV 


log 10 m L (eV) 


X 


X 


X 


-0.681 Z :lTr 


£m„ (eV) 


X 


X 


X 


o.65i 8: 8 4 2 3 9 5 


77l0/^ 777„ 


X 


X 


X 


o.337 8:i8 


mpp (eV) 


X 


X 


X 


fl °- 243 

U.lOO 0.0842 


0,/i 2 


X 


X 


X 


n lie 0.121 
U.110 0.108 


n s 


X 


X 


X 


0.949 8:118 



TABLE II: The mean value and 95% lower and upper credible intervals for various parameters and combinations of data. 



In Figs. [T] and [2] we show the likelihood contours for 
cases (a) and (b) for the assumption of normal hierar- 
chy. In both cases we have used the present uncertain- 
ties on the parameters of the mixing matrix (specified 
in Eq. [5]) which means that the non-trivial behaviour 
for small m^^mpp cannot be resolved. Given future im- 
proved constraints from reactor or long baseline experi- 
ments this region will look significantly different (see e.g. 
[47l . |48| for a thorough discussion of this point). 

1. KATRIN and GERDA only 

To complete this section we have done a parameter 
study of cases (a) and (b) using only KATRIN and 
GERDA data, excluding cosmological constraints. The 
results of this can be seen in Table HTT1 For case (a) where 
there is no detection from either experiment the combi- 
nation of cosmological data with KATRIN and GERDA 
slightly strengthens the bound on parameters, but since 
they center on the same best fit value there is no marked 



difference when cosmological data is added. 

However, this changes completely when case (b) is 
studied. Here, KATRIN and GERDA data prefer a 
higher value for ^ m„ and the combination of all three 
data sets significantly shift the allowed range for all of 
the neutrino parameters. We note that other cosmolog- 
ical parameters such as Q c h 2 and n s are not affected in 
any way when the minimal ACDM model is assumed. 
This conclusion does not hold when larger cosmological 
parameter sets are used, a point discussed in the next 
subsection. 



D. Extended cosmological models 

In order to illustrate the relation between neutrino 
experiments and cosmological parameter estimation we 
have performed the same analysis as before, but now 
adding two additional cosmological parameters to the fit: 
w, the dark energy equation of state, and a s , the run- 
ning of the scalar spectral index (giving a total of 16 
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FIG. 1: 68% and 95% contours for the Cosmo+KATRIN+GERDA data, assuming best fit values of m@ = mpp = (case a) 
Inverted hierarchy is assumed. 



Parameter 


Cosmo+KATRIN+GERDA 


KATRIN+GERDA 


Normal hierarchy 


771,3,0 = 0, 771,3,3,0 = 


log 10 m L (eV) 
E™» (eV) 

m /3/zC m » 

mpp (eV) 


-1.438l|.r 

0.168o!o79 
0.260 8;? |1 
0.0317 ° ° l 9S 


-1 476 -°- 968 

1.4/0 _ 2 go 
0.157o!o789 
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0.0293 g;^ 57 


Normal hierarchy 


771,3,0 = 0.28 eV, 771,3/3,0 = 0.18 eV 


log 10 m L (eV) 
E™* (eV) 

mgp (eV) 


-0.660 zg-.SS 

0.6748;t|I 
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-0 571 -° 493 

U.OI 1 _o.670 
n 01 7 0.969 
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TABLE III: The mean value and 95% lower and upper credible intervals for various parameters and combinations of data. 



parameters in the MCMC analysis). Particularly w is 
known to be degenerate with E m v and therefore any 
independent information on E m v from experiments is 
potentially important for dark energy physics. This par- 
ticular degeneracy has been studied quite extensively in 
recent literature. The most recent example is [44| where 
the impact of a positive KATRIN detection of mp on the 
estimation of w is discussed. 

In Fig. [3] we show the degeneracy between E m„ and 
w for our case (b), assuming normal hierarchy. The cor- 
responding numbers are shown in Table [TV] The results 



confirm previous findings, i.e. that a strongly negative 
equation of state for dark energy can be compensated by 
increasing the neutrino mass [491 , |50T | . This also means 
that the allowed region of w for case (b) is shifted to more 
negative values, in this case only marginally allowing a 
cosmological constant (the ID 95% credible interval is 
— 1.39 < w < —0.99). T he p resent result compares well 
with what is obtained in [4J] , although the assumed best 
fit values are slightly different. Note also that our treat- 
ment of cosmological data is slightly different from [44[ 
because we use the full BAO correlation function instead 
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FIG. 2: 68% and 95% contours for the Cosmo+KATRIN+GERDA data, assuming best fit values of m = 0.28 eV, rapp = 0.18 
eV (case b). Inverted hierarchy is assumed. 



Parameter 


Cosmo(+ii> + a s )+KATRIN+GERDA 


Normal hierarchy 


mp,o = 0.28 eV, m p,o = 0.18 eV 


log 10 m L (eV) 
£m„ (eV) 

mpp (eV) 
w 


0.635 _q 
0.731 Jlf 
0.330 8;li| 

0.176 0.0977 
_-, 17 -1.39 
J-- 1 ' -0.99 



TABLE IV: The mean value and 95% lower and upper credible 
intervals for case (b) with the larger cosmological parameter 
space. 



of the A functional parametrisation. For comparison, in 
Fig. [4] we show the ^ m v — w degeneracy for cosmological 
data only. In this case we find the result J2 m v < 0.56 
eV and -1.22 < w < -0.88, both at 95% CJL, a result 
completely consistent with what was found in [211 ] for the 
same model and data, but using maximisation instead of 
marginalisation. Note also that in this extended model 
the best fit x 2 increases by 2.6 compared to the case 
where only cosmological data is used (compared to 3.1 
in the smaller parameter space discussed above), i.e. in 
the extended model the inconsistency between cosmology 
and the assumed positive detection from KATRIN and 
GERDA is less pronounced. 



The exercise carried out in this subsection clearly illus- 
trates why cosmological bounds on neutrino properties 
are model dependent. Note that this degeneracy would 
be even stronger if only CMB data is considered. 
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FIG. 3: 68% and 95% contours for the Cosmo(+u> + 
a s ) +KATRIN+GERD A data in the J2 m v ~ w plane, as- 
suming best fit values of mp = 0.28 eV, nipp = 0.18 eV (case 
b). Normal hierarchy is assumed. 
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FIG. 4: 68% and 95% contours for the Cosmo(+u; + a a ) data 
in the ^ m„ — w plane. Normal hierarchy is assumed. 

IV. DISCUSSION 

A detailed neutrino parameter estimation study has 
been carried out using the Markov Chain Monte Carlo 
technique with the goal of unifying the various techniques 
for measuring the absolute neutrino mass scale. The 
MCMC technique is extremely powerful in this regard 
and allows for a very fast scanning many-dimensional 
likelihood spaces. In the concrete example here we have 
used 8 parameters describing the properties of light, ac- 
tive Majorana neutrinos, and 6 further parameters which 
specify the cosmology. 

We find that for present data the combination of 
cosmological data with the upper limit on mpp from 
Heidelberg-Moscow slightly improves the existing cosmo- 
logical bound on the sum of neutrino masses. 

More interestingly we have studied the interplay be- 
tween various future constraints from cosmology, tritium 
decay and neutrinoless double beta decay. If all probes 



come up with a negative result the addition of data sets 
does not yield any radically new information. However, 
we have also studied an example in which the upcoming 
KATRIN and GERDA experiments are both assumed 
to provide tentative evidence for neutrino mass. In this 
case the combination of all three types of data allows for 
a much stronger constraint on neutrino properties than 
otherwise allowed. 

Finally we have also studied how experimental data 
from tritium decay or neutrinoless double beta decay can 
help in cosmological parameter estimation, particularly 
concerning the dark energy equation of state. 

It should be noted that in the present analysis only 
presently available cosmological data has been used. In 
the same time frame as KATRIN and GERDA new cos- 
mological data will become available and is likely to im- 
prove the cosm ologica l neutrino mass bound significantly 
(see [5lT [52l. l53l. 154 l55l. IBlil. l57l. l58l. |59| for a non-exhaustive 
list). In the somewhat longer term cosmological con- 
straints can be potentially be pushed below 0.1 eV sen- 
sitivity to 

TTiy . At the same time neutrinoless double 
beta decay experiments will have equally improved sensi- 
tivity and it will very likely be possible to determine the 
absolute neutrino mass as well as the nature of the mass 
hierarchy. 

In conclusion, the combination of cosmological data 
with experimental neutrino data in a global analysis will 
be extremely useful in the future, when more precise ex- 
perimental data becomes available. 
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